#include "common.hpp"
#include "integration.hpp"

/** This function is designed to compute the value of f(x). Here, f is an
  * application which is consisited(by multiplication or by division) of basic functions.
  * For example:
  * 	f(x)=c*sin(x)/cos(x)*(x^5)
  * The function name is expressed by 2 int array and 1 double number.
  */
double fun_callback(function_info *funInfo, double x)
{
	/*Analyze the function's structure, and compute the function value step by step. */
	int N = funInfo->length;
	double fun_value = funInfo->coef;
	for(int i=0; i<N; i++)   
	{
		switch((funInfo->multi_fun)[i])
		{
			case SIN:
				if ((funInfo->opera_flag)[i]==MUL) 
				{   		
					fun_value = fun_value*sin(x);
					break;
				}

				else
				{
					fun_value = fun_value/sin(x);
					break;
				}
			case COS:
				if ((funInfo->opera_flag)[i]==MUL)
				{
					fun_value = fun_value*cos(x);
					break;
				}
				else
				{	
					fun_value = fun_value/cos(x);
					break;
				}
			case EXP:
				if ((funInfo->opera_flag)[i]==MUL)
				{
					fun_value = fun_value*exp(x);
					break;
				}
				else
				{	
					fun_value = fun_value/exp(x);
					break;
				}
			case LN:
				if ((funInfo->opera_flag)[i]==MUL)
				{
					fun_value = fun_value*log(x);
					break;
				}
				else
				{	
					fun_value = fun_value/log(x);
					break;
				}
			case LOG:
				if ((funInfo->opera_flag)[i]==MUL)
				{
					fun_value = fun_value*(log(x)/log(10.0));
					break;
				}
				else
				{	
					fun_value = fun_value/(log(x)/log(10.0));
					break;
				}
			default:
				if (((funInfo->opera_flag)[i])==MUL)
				{				
					fun_value = fun_value*pow(x,(funInfo->multi_fun)[i]);
					break;
				}
				else
				{ 
					fun_value = fun_value/pow(x,(funInfo->multi_fun)[i]);
					break;
				}
		}
	}
	return fun_value;
}

/** Compute the Riemann integration by accumulating the rectangular area.
  * Here, we use the left point to compute the height.
  */
double left_rect_method(function_info *funInfo,
						 double lower_limit=0.0, double upper_limit=1.0, double step=0.01)
{
	
	/* Compute the number of times for accumulation, 
	"+step/10.0" is used to avoid to getting wrong N in Windows */
	int N;
	N = (upper_limit+step/10.0-lower_limit)/step;          	
	
	double integration_val=0;			/* The integration value, initialized by 0 */
	double tmp=0;							  	
	double x=lower_limit-step;

	/* Start to accumulate the area */
	int i;
	for (i=0; i<N; i++)
	{
		x += step;

		/* Use fun_callback() to compute the function value, and accumulate the values */
		tmp += fun_callback(funInfo, x);
	}

	return step*tmp;    				/* Return the integration value */
}


/** Implementation of the trapezoid method.
  * There is strong relation between the trapezoid method and the left rectangle method,
  *	We will make use of it to simplify the code.	
  */
double trapezoid_method(function_info *funInfo,
						 double lower_limit=0.0, double upper_limit=1.0, double step=0.01)
{
	/* Get  the left_left_rect_value by left rectangle method */
	double left_rect_value=0, trapezoid_value=0;
	left_rect_value = left_rect_method(funInfo, lower_limit, upper_limit,  step);			

	/* Compute the trapezoid integration value */
	trapezoid_value = left_rect_value 
			+ (fun_callback(funInfo, upper_limit)
			- fun_callback(funInfo, lower_limit))/2*step; 		

	return trapezoid_value;
}


/** Implementation of Gauss's 5-point  method.	
  */
double gauss_5p_method(function_info *funInfo, double lower_limit=0.0, double upper_limit=1.0)
						 
{
	double cst1 = (upper_limit-lower_limit)/2.0;
	double cst2 = (upper_limit+lower_limit)/2.0;
	
	/* Get the five points */
	double point1, point2, point3, point4, point5;
	point1 = cst1*(-0.9061798)+cst2;
	point2 = cst1*(-0.5384693)+cst2;
	point3 = cst1*(-0.0000000)+cst2;
	point4 = cst1*(+0.5384693)+cst2;
	point5 = cst1*(+0.9061798)+cst2;

	/* Use the ready-made formula to compute the integration value */
	double gauss_5p_value = cst1*( 0.2369269*fun_callback(funInfo, point1)
								   +0.4786287*fun_callback(funInfo, point2)
								   +0.5688889*fun_callback(funInfo, point3)
								   +0.4786287*fun_callback(funInfo, point4)
								   +0.2369269*fun_callback(funInfo, point5));
	return gauss_5p_value;
}

/** Implementation of the Monte Carlo method.	
  */
double mc_method(function_info *funInfo, double lower_limit=0.0, double upper_limit=1.0,
						int num_of_experiment=1000, int num_of_part=1000)
{
	double part;
	part = (upper_limit-lower_limit)/num_of_part;  			/* Generate a partition */
	double integration_val = 0;					/* The integration value */


	for(int i=1; i<=num_of_part; i++)
	{
		
		int positive_counter = 0;				/* Record the number that random point drops in the positive area */
		int negative_counter = 0;				/* Record the number that random point drops in the negative area */
		int invalid_exp_counter = 0;			/* Record the number of invalid experiments */


		/* Look for a appropriate rectangle */
		double left_point = lower_limit+(i-1)*part;
		double right_point = lower_limit+i*part;
		double extrenum[2]={0,0};
		get_ymin_ymax(funInfo, left_point, right_point, extrenum, 0.001);



		for(int j=1; j<=num_of_experiment; j++)
		{
			/* Generate a 2D random point */
			double rand_x = random(left_point, right_point);
			double rand_y = random(extrenum[0], extrenum[1]);
			
			/* Random point drops in the positive area */
			if(rand_y>0 && rand_y<fun_callback(funInfo, rand_x) )  	
				positive_counter++;	
			
			/* Random point drops in the negative area */												
			if(rand_y<0 && rand_y>fun_callback(funInfo, rand_x))   	
				negative_counter++;

			/* Random point drops in  the border line */
			if(rand_y=0 || rand_y==fun_callback(funInfo, rand_x)) 	
						invalid_exp_counter++;
		}

		/* Compute the frequence */
		double freq = ((double) positive_counter-negative_counter)/((double) num_of_experiment-invalid_exp_counter); 

		/* Compute the integration*/
		integration_val += freq*part*(extrenum[1] - extrenum[0]);	
	}

	return integration_val;
}


/* A fuction for finding an appropriate weight of rectangle */
void get_ymin_ymax(function_info *funInfo, double lower_limit, double upper_limit, double extrenum[], double step=0.01)
{
	int N = (upper_limit - lower_limit)/step;
	for(int i=0; i<=N; i++)
	{
		double x = lower_limit+i*step;
		double f_x = fun_callback(funInfo, x);

		if (f_x <= extrenum[0])
			extrenum[0] = f_x;
		if (f_x >= extrenum[1])
			extrenum[1] = f_x;
	}

	/** What we have got are approximate extrenums,
	  * so we enlarge them to make sure they are valid */
	extrenum[0] = 1.2*extrenum[0];
	extrenum[1] = 1.2*extrenum[1];
}


/* Generate a random number in a given interval */
double random(double lower_limit, double upper_limit)
{
	double rand_num;
	rand_num = (double) (rand());
	rand_num = rand_num/RAND_MAX;
	rand_num = lower_limit + (upper_limit - lower_limit)*rand_num;
	return rand_num;
}
